{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "mSIE7vZ2hpKJ"
      },
      "source": [
        "Licensed under the Apache License, Version 2.0"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Ewkfmus07yMq"
      },
      "outputs": [],
      "source": [
        "#@title Install dependencies\n",
        "!pip install xarray\n",
        "!pip install gcsfs\n",
        "!pip install pyresample"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "_Ry4jMdu7z7k"
      },
      "outputs": [],
      "source": [
        "import datetime\n",
        "import pprint\n",
        "\n",
        "import gcsfs\n",
        "import matplotlib.pyplot as plt\n",
        "import numpy as np\n",
        "import pyresample\n",
        "import xarray"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "SCzJYBRqsmxs"
      },
      "source": [
        "# Load GOES data\n",
        "\n",
        "GOES filenames contain lots of information, with a few important parts.\n",
        "\n",
        "`OR_ABI-L1b-RadF-M6C11_G16_s20230010020206_e20230010029514_c20230010029565.nc`\n",
        "* `L1b` - Level-1b data product, as opposed to L2 which are derived data.\n",
        "* `RadF` - Full-disk radiances.\n",
        "* `C11` - Data for band #11.\n",
        "* `s20230010020` - Start time of the GOES scan. Scans happen every 10 minutes, so this scan happened on:\n",
        "  * `2023` - year\n",
        "  * `001` - day of year, e.g. Jan 1st\n",
        "  * `00` - hour, `20` - minute\n",
        "\n",
        "Here we load 3 bands from this scan - `11`, `14`, and `15` and convert the radiances to brightness temperatures. These 3 bands are required for our color scheme.\n",
        "\n",
        "GOES data is hosted on Google Cloud for public use.\n",
        "\n",
        "References:\n",
        "* GOES data products: https://www.goes-r.gov/products/overview.html\n",
        "* GOES public dataset: https://console.cloud.google.com/marketplace/product/noaa-public/goes-16"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "ZPuI9M2_8Fda"
      },
      "outputs": [],
      "source": [
        "fs = gcsfs.GCSFileSystem(project='gcp-public-data-goes-16')\n",
        "print('All band paths for scan:')\n",
        "pprint.pprint(fs.glob('gcp-public-data-goes-16/ABI-L1b-RadF/2023/001/00/*_s20230010020*'))\n",
        "\n",
        "paths = {\n",
        "    11: 'gcp-public-data-goes-16/ABI-L1b-RadF/2023/001/00/OR_ABI-L1b-RadF-M6C11_G16_s20230010020206_e20230010029514_c20230010029565.nc',\n",
        "    14: 'gcp-public-data-goes-16/ABI-L1b-RadF/2023/001/00/OR_ABI-L1b-RadF-M6C14_G16_s20230010020206_e20230010029514_c20230010029555.nc',\n",
        "    15: 'gcp-public-data-goes-16/ABI-L1b-RadF/2023/001/00/OR_ABI-L1b-RadF-M6C15_G16_s20230010020206_e20230010029522_c20230010029573.nc'\n",
        "}\n",
        "\n",
        "brightness_temperatures = {}\n",
        "for band_id, path in paths.items():\n",
        "  with fs.open(path, 'rb') as f:\n",
        "    dataset = xarray.open_dataset(f)\n",
        "    # Convert radiances to brightness temperature\n",
        "    brightness_temperatures[band_id] = (dataset.planck_fk2.data / np.log((dataset.planck_fk1.data / dataset.Rad.data) + 1) - dataset.planck_bc1.data) / dataset.planck_bc2.data\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bdCFygO-wbP-"
      },
      "source": [
        "# Generate false color image\n",
        "\n",
        "In order to view contrails in GOES, we use the \"ash\" color scheme. This color scheme was originally developed for viewing volcanic ash in the atmosphere but is also useful for viewing thin cirrus, including contrails. In this color scheme, contrails appear in the image as dark blue.\n",
        "\n",
        "Note that we use a modified version of the ash color scheme here, developed by Kulik et al., which uses slightly different bands and bounds tuned for contrails.\n",
        "\n",
        "References:\n",
        "* Ash Color Scheme (page 7): https://eumetrain.org/sites/default/files/2020-05/RGB_recipes.pdf"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "gH62mImh-Skd"
      },
      "outputs": [],
      "source": [
        "_T11_BOUNDS = (243, 303)\n",
        "_CLOUD_TOP_TDIFF_BOUNDS = (-4, 5)\n",
        "_TDIFF_BOUNDS = (-4, 2)\n",
        "\n",
        "def normalize_range(data, bounds):\n",
        "  \"\"\"Maps data to the range [0, 1].\"\"\"\n",
        "  return (data - bounds[0]) / (bounds[1] - bounds[0])\n",
        "\n",
        "def false_color_image(brightness_temperatures):\n",
        "  \"\"\"Generates ash false color image from GOES brightness temperatures.\"\"\"\n",
        "  r = normalize_range(brightness_temperatures[15] - brightness_temperatures[14], _TDIFF_BOUNDS)\n",
        "  g = normalize_range(brightness_temperatures[14] - brightness_temperatures[11], _CLOUD_TOP_TDIFF_BOUNDS)\n",
        "  b = normalize_range(brightness_temperatures[14], _T11_BOUNDS)\n",
        "  return np.clip(np.stack([r, g, b], axis=-1), 0, 1)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "4RUMqNr7AF7j"
      },
      "outputs": [],
      "source": [
        "fci = false_color_image(brightness_temperatures)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "c4Fxbb6EFLrN"
      },
      "outputs": [],
      "source": [
        "plt.figure(figsize=(10, 10))\n",
        "plt.imshow(fci)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "s-oQINSFCxdU"
      },
      "source": [
        "# Find a location on the image\n",
        "\n",
        "The GOES data comes with information about the projection, which we can use to find a lat/lng location on the image. Note that we are *not* considering parallax here, so the lat/lng location that we plot will be on the surface.\n",
        "\n",
        "In this example, we find New York City on the image."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "HAeQB40AJkqL"
      },
      "outputs": [],
      "source": [
        "NYC_LAT = 40.7128\n",
        "NYC_LNG = -74.0060"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "qz6lqn1KFMdv"
      },
      "outputs": [],
      "source": [
        "# Open the netCDF file for one of the bands. All bands will have the same projection information.\n",
        "with fs.open(paths[11], 'rb') as f:\n",
        "  dataset = xarray.open_dataset(f)\n",
        "\n",
        "  h0 = dataset.goes_imager_projection.perspective_point_height\n",
        "  area_def = pyresample.geometry.AreaDefinition(\n",
        "    area_id='all_goes_16',  # Used only for pyresample logging\n",
        "    proj_id='deprecated',  # Deprecated but required by pyresample\n",
        "    description='all_goes_16',  # Used only for pyresample logging\n",
        "    projection={  # proj4 dict\n",
        "        'proj': 'geos',  # Stands for 'geostationary'\n",
        "        'units': 'm',\n",
        "        'h': str(h0),\n",
        "        'lon_0': str(\n",
        "            dataset.goes_imager_projection.longitude_of_projection_origin\n",
        "        ),\n",
        "        'a': str(dataset.goes_imager_projection.semi_major_axis),\n",
        "        'b': str(dataset.goes_imager_projection.semi_minor_axis),\n",
        "        'sweep': dataset.goes_imager_projection.sweep_angle_axis,\n",
        "    },\n",
        "    width=dataset['x'].shape[0],\n",
        "    height=dataset['y'].shape[0],\n",
        "    area_extent=[\n",
        "        dataset['x_image_bounds'].data[0] * h0,\n",
        "        dataset['y_image_bounds'].data[1] * h0,\n",
        "        dataset['x_image_bounds'].data[1] * h0,\n",
        "        dataset['y_image_bounds'].data[0] * h0,\n",
        "    ],\n",
        "  )\n",
        "# Convert the lat/lng to row/col indices in the image.\n",
        "col, row = area_def.lonlat2colrow(NYC_LNG, NYC_LAT)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "4y3umWwrJVoS"
      },
      "outputs": [],
      "source": [
        "plt.figure(figsize=(10, 10))\n",
        "plt.plot(col, row, marker='x', color='red', markersize=12, markeredgewidth=5)\n",
        "plt.imshow(fci)\n",
        "plt.show()"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "private_outputs": true,
      "provenance": []
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}
